

n =64;
N = n*n;

%samplePIC = SampPIC(20,20,20,20,n);



%set projection num
 L = 180;
 
%measure matrix
matrix = getSampMatrix(n,L);
A = double(matrix);

%static sample
% b = A * samplePIC(:);
% bback = b;

theta = 1:(180/L):180;
%[R] = radon(samplePIC,theta);
%Rback = R;
%xp = iradon(R,theta,n);
%xpStatic = xp(:);

%dyna sample
cycleNum = 5;
cycleLen = floor(L/cycleNum);
halfCycleLen = floor(cycleLen/2);
MaxDis = 4;
speed = MaxDis/halfCycleLen;
ObjSize = 4;
if halfCycleLen == 0
    disp('need more projection')
    return;
end

for i = 1:L
    j = mod(i,cycleLen);
    pos = j - halfCycleLen;
    dis = pos * speed;
    dis = floor(dis);
    a(j+1) = dis;
   
    samplePIC = SampPIC(30+dis,30+ObjSize+dis,30,30+ObjSize,n);

    Rtemp = radon(samplePIC,theta(i));
    R(:,i) = Rtemp;
    Ttheta(j+1,ceil(i/cycleLen)) = i;
    
    btemp = A * samplePIC(:);
   %-----------------add noise----
%    noise =   randn(size(btemp));
%    btemp = btemp +    25 * noise/norm(noise);
	%-------------------------------
        proj(j+1,(ceil(i/cycleLen)-1)*n+1:(ceil(i/cycleLen))*n) = (i-1)*n+1:i*n;
        bcycle(j+1,(ceil(i/cycleLen)-1)*n+1:(ceil(i/cycleLen))*n) = btemp((i-1)*n+1:i*n);
    %b((i-1)*n+1:i*n,:)= btemp((i-1)*n+1:i*n,:);       
end

xpDyna = iradon(R,theta,'spline','Shepp-Logan',n);




%----------reconstruction-------------------------------------
%set lamda
lamda = 0;

Phase = [1:2];
%
detVal =   bcycle(Phase,:)';
detVal = detVal(:);
PrjSeq = proj(Phase,:)';
PrjSeq = PrjSeq(:);
thetaSeq = Ttheta(Phase,:)';
thetaSeq = thetaSeq(:);

cvx_begin
    variable xrec(N)
    %minimize( norm(xrec,1) )  
    minimize(lamda*norm(xrec - xpDyna(:),1)+ (1-lamda)*norm(xrec,1))
    %minimize (lamda*sum_square(xrec-xpDyna(:)) + (1-lamda)*norm(xrec,1))
    subject to
     %   bcycle(5,0*n+1:5*n)'  == A(proj(5,0*n+1:5*n),:) * xrec ;   
     %detVal  == A(PrjSeq,:) * xrec ;   
     norm(detVal  - A(PrjSeq,:) * xrec) < 25
cvx_end

%R = radon(samplePIC,theta(Ttheta(10,:)));
bp= iradon(R(:,thetaSeq),theta(thetaSeq),64);

figure(1)
subplot(3,1,1)
showPIC(xpDyna(:),n,n);
subplot(3,1,2);
showPIC(bp(:),n,n);
subplot(3,1,3)
showPIC(xrec,n,n);

%{
 figure(1)
subplot(2,2,1)
showPIC(xpDyna(:),n,n);
title('Prior Image');
subplot(2,2,2);
showPIC(bp(:),n,n);
title('backprojection reconstructed Image');
subplot(2,2,3)
showPIC(xrec,n,n);
title('CS reconstructed Image');
subplot(2,2,4)
showPIC(xback,n,n);
title('PICCS reconstructed Image');
%}